Hyperlink to webpage:

https://Carolinah23.github.io

Final Tutorial - Analysing data from Water Wells on Baton Rouge, LA

Diana Carolina Hurtado Pulido
December/16/2021

CMPS 3160/6160

Description project and Goals

Currently, I am working on finding subsidence rates in Baton Rouge, Louisiana during the last two decades, and investigating what factors are causing subsidence using LiDAR (Light Detection And Ranging) data from 1999 and 2018.

Subsidence is the vertical land movement caused by natural and anthropogenic factors such as sediment compaction, isostatic adjustments, fault slip, and extraction or injection of fluids. The Gulf of Mexico coastline is under constant monitoring due to high rates of Sea Level rise and subsidence, which cause rapid land loss. The study area is not in the coastal area though it is subsiding (Figure 1). This area is of particular interest because there has not been enough research to determine what factors are causing subsidence. This area has two geological faults, grow of urban development, almost 2000 water wells (with different uses) active during the study period, approximately 40 active oil and gas wells, and 11 injection wells.

So far, my results show subsidence in the whole region (Figure 1). Surprisingly subsidence increases from south to north, which is the opposite of the expected results. Areas closer to the coast and water bodies have large subsidence rates due to younger sediments compaction, but most important in this area, the faults move towards the south. Interestingly, small areas show localized subsidence and uplifting, this behavior likely is related to human activities. Then, having these results, the main goal of this project is to find out how these subsidence values are related to groundwater extraction. Important questions:

1. The location of clusters of wells effectively occurs in areas subsiding or uplifting locally?
2. Are deeper wells causing more or subsidence than wells at shallow depths?
3. Is there any particular well use that may be causing more vertical changes?

Figure 1: Relative subsidence on Baton Rouge using LiDAR differencing between 1999 and 2018. Each map shows different methods applied to find elevation changes (Result of my research)

from IPython.display import Image
Image(filename="../../Images/GCD_raws.jpg")

Extraction, Transform and Load (ETL)

Data description

Dataset 1 and 2 come from the Department of Natural Resources of Louisiana. The dataset 1 - Wells_df - correspond to the water wells functioning during the period of study (1999 and 2018) in the area. We collected this data (with a collaborator), last year (2020), and calculated the well depth on meters. The data is not complete, the missing values were not published when the well-owners uploaded the information, or the information is too old in some cases. Dataset 2 locations_df are the coordinates on NAD83(2011) / UTM zone 15N and other data from the water wells, just some general information will be conserved for this analysis.
Dataset 3 - varZ_df has the information of elevetation changes calculated in my research between 1999 and 2018 in a grid of 100 meters in the same coordinate system than locations_df .

For my analysis I will use the following variables in Wells_df: Well Depth (on meters), Well Use, Yield (rate of extraction - gallons per minute), and the dates of construction, last date the well was active..
For the dataset locations_df: Coordinates X and Y, LocalWellNumber (to merge the data with the first dataset), Water table depth (how depth can we find underground water - in meters), and Aquifer Name.
For dataset varZ_df : Coordinates X and Y, and average elevation change in each point.

Loading data:

I and uploading the data using the read_csv()

#Importing libraries
from bs4 import BeautifulSoup
import re
import requests as rq
import pandas as pd
import numpy as np
import lxml
import datetime as dt
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.neighbors import KNeighborsRegressor
from sklearn.feature_extraction import DictVectorizer
from sklearn.preprocessing import StandardScaler

#Loading data
Wells_df=pd.read_csv("../data/Water_Well_Data.csv")
locations_df=pd.read_csv("../data/Locations.csv")
varZ_df=pd.read_csv("../data/Subsidence_19years_Ground.csv")

Wells_df.head()

	LocalWellNumber	WellDepth	Well_Depth_Meters	WellUse	Active	LastActive_Plugged	LastActive_Plugged_date	DateConstructed	Yield	SerialNumber
0	271	650	198.12	Domestic	Y	Current	08/02/21	1/1/1919	NN	NN
1	368	450	137.16	Irrigation	Y	Current	08/02/21	10/1/1992	30	NN
2	621	1487	453.2376	Public Supply	Y	Current	08/02/21	12/1/1956	295	1033005
3	804A	1950	594.36	Observation	Y	Current	08/02/21	2/1/1966	30	NN
4	804B	2762	841.8576	Observation	Y	Current	08/02/21	2/1/1966	25	NN

locations_df.head()

	LOCAL_WELL	OWNERS_NAM	DATE_COMPL	GEOLOGIC_U	WATER_LEVE	DATE_MEASU	PUMP_SETTI	PUMP_RATE	YIELD_MEAS	PUMP_HOURS	PUMP_DOWN_	AQUIFER_NA	X	Y
0	10180Z	KEAN'S THE CLEANERS	10/16/2020	112SESC	NaN	NaN	NaN	NaN	NaN	NaN	G	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	680528	3361312
1	9749Z	BFI ENERGY SYSTEMS OF BOSTON, INC	5/7/12	112SESC	13.0	5-7-12	NaN	NaN	NaN	NaN	P	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	686884	3364377
2	9750Z	BFI ENERGY SYSTEMS OF BOSTON, INC	5/7/12	112SESC	8.0	5-7-12	NaN	NaN	NaN	NaN	P	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	686883	3364408
3	9751Z	BFI ENERGY SYSTEMS OF BOSTON, INC	5/7/12	112SESC	6.0	5-7-12	NaN	NaN	NaN	NaN	P	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	686830	3364407
4	9752Z	BFI ENERGY SYSTEMS OF BOSTON, INC	5/7/12	112SESC	6.0	5-7-12	NaN	NaN	NaN	NaN	P	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	686857	3364407

varZ_df.head()

	fid	VarZ	X	Y
0	1	0.0	675085	3376455
1	2	0.0	675185	3376455
2	3	0.0	675285	3376455
3	4	0.0	675385	3376455
4	5	0.0	675485	3376455

Tidy data

1. Droping
2. Raplacing incorrect or empty values by NaN
3. Convert data to appropiate types
4. Renaming and organizing columns

Wells_df

Wells_df has two variables that I will not use: WellDepth and SerialNumber. The first variable will not be used because I will use the information on meters, and the second one is not the main identificator of the observations and is not complete.

Wells_df=Wells_df.drop(["WellDepth","SerialNumber", "LastActive_Plugged", "Yield"], axis=1)

There are "NN" values in the LastActive_Plugged_date and DateConstructed that must be changed to NaN using np.nan

Wells_df["LastActive_Plugged_date"].replace("NN", np.nan, inplace=True)
Wells_df["DateConstructed"].replace("NN", np.nan, inplace=True)

The cell below shows the variable types, to datetime the variables Last_Active_Plugged_date and DateConstructed and to number I changed the variable Well_Depth_Meters

Wells_df.dtypes

LocalWellNumber            object
Well_Depth_Meters          object
WellUse                    object
Active                     object
LastActive_Plugged_date    object
DateConstructed            object
dtype: object

Wells_df["Well_Depth_Meters"]=pd.to_numeric(Wells_df["Well_Depth_Meters"], errors="coerce")

Wells_df["LastActive_Plugged_date"]=pd.to_datetime(Wells_df["LastActive_Plugged_date"], format="%m/%d/%y")
Wells_df["DateConstructed"]=pd.to_datetime(Wells_df["DateConstructed"], format="%m/%d/%Y")

For future models, is better to have the time that a well was active. This variable is the differente between LastActive_Plugged_date and DateConstructed variables.

Wells_df["ActiveDays"]=(Wells_df["LastActive_Plugged_date"]-Wells_df["DateConstructed"]).dt.days

Wells_df.dtypes

LocalWellNumber                    object
Well_Depth_Meters                 float64
WellUse                            object
Active                             object
LastActive_Plugged_date    datetime64[ns]
DateConstructed            datetime64[ns]
ActiveDays                        float64
dtype: object

And finally, I am replacing the Y and N (Yes and No) by 1 and 0 in the Active column. This variable indicates if the well was active extracting groundwater to the date that we took the data. For instance, if the well was active between 2001 and 2014 it should has a N, also, if the well is been used to monitior groundwater it should be inactive.

Wells_df["Active"].replace("Y", 1, inplace=True)
Wells_df["Active"].replace("N", 0, inplace=True)

Wells_df.head()

	LocalWellNumber	Well_Depth_Meters	WellUse	Active	LastActive_Plugged_date	DateConstructed	ActiveDays
0	271	198.1200	Domestic	1	2021-08-02	1919-01-01	37469.0
1	368	137.1600	Irrigation	1	2021-08-02	1992-10-01	10532.0
2	621	453.2376	Public Supply	1	2021-08-02	1956-12-01	23620.0
3	804A	594.3600	Observation	1	2021-08-02	1966-02-01	20271.0
4	804B	841.8576	Observation	1	2021-08-02	1966-02-01	20271.0

locations_df

locations_dfFirst, I will calculate the water table depth o meters because it is on feet. Having this variable I will drop a set of variables that are not necessary in this analysis to conserve only the variables mentioned above.

locations_df["Water_depth"]=(locations_df["WATER_LEVE"]/0.3048)
locations_df=locations_df.drop(["OWNERS_NAM", "DATE_COMPL", "GEOLOGIC_U", "DATE_MEASU", "PUMP_SETTI", "PUMP_RATE", 
                                "YIELD_MEAS", "PUMP_HOURS", "PUMP_DOWN_", "WATER_LEVE"], axis=1)

Now, I change the name of the identificator with rename to have the same column in both tables, also I reorganize the columns for better understanding.

locations_df=locations_df.rename({"LOCAL_WELL":"LocalWellNumber"}, axis=1)
locations_df=locations_df[["X", "Y", "LocalWellNumber", "Water_depth", "AQUIFER_NA"]]

locations_df.head()

	X	Y	LocalWellNumber	Water_depth	AQUIFER_NA
0	680528	3361312	10180Z	NaN	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...
1	686884	3364377	9749Z	42.650919	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...
2	686883	3364408	9750Z	26.246719	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...
3	686830	3364407	9751Z	19.685039	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...
4	686857	3364407	9752Z	19.685039	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...

varZ_df

varZ_df=varZ_df.drop("fid", axis=1)
varZ_df=varZ_df[["X", "Y", "VarZ"]]

varZ_df.head()

	X	Y	VarZ
0	675085	3376455	0.0
1	675185	3376455	0.0
2	675285	3376455	0.0
3	675385	3376455	0.0
4	675485	3376455	0.0

Exploratory Data Analysis

For this second milestone, I am doing some analysis to see how my data distributes using different variables.

First I will merge the tables with the information about the wells with its respective localization. These dataframes have different size because locations_df has the information of all the wells that have been in the area since the data started to be collected, while Wells_df only has data for wells that were active in any period between 1999 and 2018.

Wells_located_df=locations_df.merge(Wells_df, on=["LocalWellNumber"])
Wells_located_df.head()

	X	Y	LocalWellNumber	Water_depth	AQUIFER_NA	Well_Depth_Meters	WellUse	Active	LastActive_Plugged_date	DateConstructed	ActiveDays
0	680528	3361312	10180Z	NaN	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	4.5720	Monitor	1	2021-08-02	2020-10-16	290.0
1	686884	3364377	9749Z	42.650919	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	9.1440	Recovery	1	2021-08-02	2012-05-07	3374.0
2	686883	3364408	9750Z	26.246719	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	8.8392	Recovery	1	2021-08-02	2012-05-07	3374.0
3	686830	3364407	9751Z	19.685039	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	9.1440	Piezometer	1	2021-08-02	2012-05-07	3374.0
4	686857	3364407	9752Z	19.685039	SOUTHEAST LOUISIANA AQUIFER SYSTEM SURFICIAL C...	9.1440	Piezometer	1	2021-08-02	2012-05-07	3374.0

Now, I want to know how is the distribution of wells depth and water table depth. My first guess is that they should follow a similar distribution, because if the wells have the purpose of extracting water, then they should be at least as depth as the water table or maybe a little more.
The data of both variables are skewed to the right, then I applied a transformation to the variables in the graph using the Ladder of Powers with an exponent of 0.3 in both cases, then if a well has a depth of 10 meters it will appear in the graph as 1.995...

fig, ax=plt.subplots(1,2, figsize=(10,5), sharey=True)
fig.suptitle("Fig 2 - Water well depths vs Water table top")
ax[0].hist((Wells_located_df["Well_Depth_Meters"]**0.3), color="#985d00", alpha=0.7, bins=6, label="Well depth")
ax[0].legend(loc="upper right")
ax[1].hist((Wells_located_df["Water_depth"]**0.3), color="#0086ff", alpha=0.7, bins=6, label="Water table top")
ax[1].legend(loc="upper right")
fig.text(0.46, 0.05, "Meters^0.3")
plt.show()

The graph above shows that both variables do not have distributions as similar as I thought. Water table depth accumulates between 0 and 4 (0-102 meters depth), but the majority of wells just reach 40 meters, meaning that many do not reach the water table. Maybe the wells that do not reach the water table are used for other purposes, or there is always the possibility that one (or both) datasets were not filled initially correctly. We found that data for wells is not well stored. Definitely, there are more water wells at shallower depths where the sediment may be less compacted. The study area is on stable sediments, but changes in volume are still possible.

---

Now, I want to know how the wells' installation and deactivation are distributed over time, are old wells still active? when were the wells installed and how much have they been there?
The following graph shows that the construction of wells increased greatly in the '80s and decreased recently, however many of these wells are still active. Comparing the wells installed before 1980 and the wells the last active between 1980 and 1990, I can see that the rate of installation and inactivation is not even close.

#instalation vs last active date
fig, ax=plt.subplots(figsize=(9,5))
fig.suptitle("Fig 3 -Date construction and last active date for wells in the area")
plt.hist(Wells_located_df["DateConstructed"], color="#ce7e00", alpha=0.7, bins= 10,label="Well construction date")
plt.hist(Wells_located_df["LastActive_Plugged_date"], color="#060757", alpha=0.7, bins= 5, label="Last date where the wells were active")
plt.legend(loc="upper left")
plt.show()

Time_Active=(Wells_located_df["LastActive_Plugged_date"]-Wells_located_df["DateConstructed"])
print("Oldest active well has been there during", Time_Active.max())
print("Average days that a well is active", Time_Active.mean())

Oldest active well has been there during 38505 days 00:00:00
Average days that a well is active 7254 days 19:08:17.423887488

---

The following graph shows how depth are the wells that exist (or have been active during the study period) for each use.
Public supply has the deepest wells, these wells provide water to water and sewage users in East Baton Rouge. The next ones are Observation wells and Industrial Wells. Then most of the wells in these cases are on areas where the water table is deep even if they are uncommon (first graph).

Wells_Uses_depth=Wells_located_df.groupby("WellUse").Well_Depth_Meters.mean()
#Well use vs mean well depth
fig=plt.figure(figsize=(10,5))
fig.suptitle("Fig 4 - Whats the average depth for each type of well use?")
Wells_Uses_depth.plot.bar(color="#386c22", edgecolor="black", width=0.8)
plt.ylabel("Mean depth (m)")
plt.show()

---

Now, I want to see how subsidence change from south to north and from west to east. This graph will help us to define what are the spatial trends that may not be easily observable in the map shown above.

The graph shows that there is more variation in elevation changes from south to north than from west to east. Just seeing the top graph, the northern area reaches variations of -0.25 m (25 centimeters), which is unusual having that the southern area should be having more negative elevations changes. With respect to the horizontal coordinate, elevation variation seems to be more stable, except for the strong bumps at the coordinate 687500, meaning that we have more variation of subsidence to the south than horizontally. This conclusion has also been studied with GPS permanent stations in the area.

VarZ_Y=varZ_df.groupby("Y").VarZ.mean() #Mean elevation change on Y
VarZ_X=varZ_df.groupby("X").VarZ.mean() #Mean elevation change on X

#Profile mean values from north to south
fig, ax=plt.subplots(2,1, figsize=(10,10))
fig.suptitle("Fig 5 - Average elevation change between 1999 and 2018")
ax[0].plot(VarZ_Y.index, VarZ_Y, color="#e69138")
ax[0].set_ylabel("Mean elevation change (m)")
ax[0].set_xlabel("Coordinate on Y (m)")
ax[0].text(3358000, -0.3, "South", backgroundcolor="#bcbcbc")
ax[0].text(3376000, -0.3, "North", backgroundcolor="#bcbcbc")
ax[1].plot(VarZ_X.index, VarZ_X)
ax[1].set_ylabel("Mean elevation change (m)")
ax[1].set_xlabel("Coordinate on X (m)")
ax[1].text(675000, -0.45, "West", backgroundcolor="#bcbcbc")
ax[1].text(694000, -0.45, "East", backgroundcolor="#bcbcbc")
plt.show()

---

In the following graph, I want to see how is the distribution of the wells in the area using the coordinates X and Y, and how depth they are. In the top graph, we see that deeper wells are located to the north, where there is more subsidence (strong negative elevation changes - graph above). Most of the wells have a depth of less than 200 meters, then, is this a hint pointing that deeper wells relate to higher values of subsidence? There are more wells in the northern area than the south, which may also indicate this.
On the other hand, the location of the wells in the X coordinate does not seem to indicate any relation with elevation changes previously seen.

#Profile mean values from north to south
fig, ax=plt.subplots(2,1, figsize=(10,10))
fig.suptitle("Fig 6 - Wells distribution on area and depth")
ax[0].scatter(Wells_located_df["Y"], Wells_located_df["Well_Depth_Meters"], color="#e69138", s=15)
ax[0].set_ylabel("Well depth (m)")
ax[0].set_xlabel("Coordinate on Y (m)")
ax[0].text(3358000, -100, "South", backgroundcolor="#bcbcbc")
ax[0].text(3377500, -100, "North", backgroundcolor="#bcbcbc")
ax[1].scatter(Wells_located_df["X"], Wells_located_df["Well_Depth_Meters"], s=15)
ax[1].set_ylabel("Well depth (m)")
ax[1].set_xlabel("Coordinate on X (m)")
ax[1].text(673000, -100, "West", backgroundcolor="#bcbcbc")
ax[1].text(696000, -100, "East", backgroundcolor="#bcbcbc")
plt.show()

Model Questions

The analysis made before brings up two important relations that will be settled using statistical modeling, 1) Relationship between calculated elevation changes and location of water wells from south-north, apparently, wells at deeper depths could cause more subsidence and 2) Public supply, Industry, and observational wells are the deepest wells, may one of these uses causing more elevation changes.

To answer these questions I will do a stadistical analysis to find the correlation between depth and elevation change in the study period. I compare these results with the distribution of uses to try to point out what use is causes more changes.

Also, I will the Yield variable to find out if the extraction rate is linked in some way to subsidence or if the depth is more important. This correlation analysis is limited because only 185 wells of 1972 in the area have this information.

Modeling

For this model I will find if there is any correlation between elevation change values calculated in my research and location from North to South, and with well depth.

Secondly, I will do a regression model where I can incorporate variables as well-use and time that the well has been working.

Note: For the sake of this project I am just using well variables to explain elevation changes, but subsidence may be caused by other factors such as fault slip, compaction of sediments, isostatic adjustments, etc... and the interaction of these.
To facilitate calculations I am assigning a subsidence value to each well using its coordinates to find the closest subsidence pixel.

Variations=[]
k=0
for i in range(len(Wells_located_df)): #1885 wells
    dist=999999
    mini=0
    for j in range(len(varZ_df)): #More or less 35000 subsidence values
        newdist=1
        newdist=np.sqrt(((Wells_located_df.at[i, "X"]-varZ_df.at[j, "X"])**2)+
                        ((Wells_located_df.at[i, "Y"]-varZ_df.at[j, "Y"])**2))
        if newdist<dist:
            dist=newdist
            mini=varZ_df.at[j, "VarZ"]
    k=k+1
    Variations.append(mini)

Wells_located_df["VarZ"]=Variations

1. Correlation analysis

The results below show that there is not a strong correlation between location and well-depth and elevation changes, actually these results suggest that there is a negative correlation between the X coordinate and elevation changes. From figure 6 we can see that just some wells have large depths (>100 meters), exactly 283 wells out of 1886 (15%), as most of the wells are at shallow depths and subsidence is still happening due to different factors, these few wells are not causing a strong effect in the general picture.

Finally, as an extra metric I calculate the correlation between wells' depth and water table to have a numerical description of figure 2, and seems that these data are well correlated tough they should be more similar in my opinion.

variables=["X", "Y", "Well_Depth_Meters", "VarZ", "Water_depth", "ActiveDays"]
display(Wells_located_df[variables].corr())

#Correlation elevation change and location of wells
Corr_ZY=Wells_located_df["VarZ"].corr(Wells_located_df["Y"])

#Correlation elevation change and Well depts
Corr_ZDepth=Wells_located_df["VarZ"].corr(Wells_located_df["Well_Depth_Meters"])

print("Correlation Elevation change and north coordinate of water wells", Corr_ZY)
print("Correlation Elevation change and water wells depth", Corr_ZDepth)

	X	Y	Well_Depth_Meters	VarZ	Water_depth	ActiveDays
X	1.000000	-0.200829	0.023632	-0.130676	-0.078027	0.001067
Y	-0.200829	1.000000	0.051364	0.031299	0.126591	0.086578
Well_Depth_Meters	0.023632	0.051364	1.000000	0.068949	0.723933	0.537543
VarZ	-0.130676	0.031299	0.068949	1.000000	0.067625	0.046391
Water_depth	-0.078027	0.126591	0.723933	0.067625	1.000000	0.298482
ActiveDays	0.001067	0.086578	0.537543	0.046391	0.298482	1.000000

Correlation Elevation change and north coordinate of water wells 0.03129919613707609
Correlation Elevation change and water wells depth 0.068948724650082

Corr_DepthWater=Wells_located_df["Well_Depth_Meters"].corr(Wells_located_df["Water_depth"])
print("Correlation well depth and water table depth", Corr_DepthWater)

Correlation well depth and water table depth 0.723933262405864

len(Wells_located_df[Wells_located_df["Well_Depth_Meters"]>100])

283

2. Regression Models

For this regression model I am using almost all variables in the dataframe to describe Elevation changes VarZ. I divided the dataset Wells_located_df in two datasets of equal size and dropped the rows with nan values. I dropped those because they cause problems when runing the model. And finally, I calculated the root mean square error to measure how far are the predictions from the real data. The error is close to 0.25 (m), considering that values of elevation change should be very small I consider that this error is pretty large. From the correlation analysis I knew that the independent variables are not great to explain this phenomena.
I tried two more models with less variables.

#Training and validation data
Training=Wells_located_df.sample(frac=0.5)
Training.dropna(inplace=True)
Validation=Wells_located_df.drop(Training.index)
Validation.dropna(inplace=True)

variables=["X", "Y", "Water_depth", "AQUIFER_NA", "Well_Depth_Meters", "WellUse", "ActiveDays"]
Independent_Train_dict=Training[variables].to_dict(orient="records")
Independent_Val_dict=Validation[variables].to_dict(orient="records")

Dependent_Train=Training["VarZ"]
Dependent_Val=Validation["VarZ"]

#Dummy encoding
vec = DictVectorizer(sparse=False)
vec.fit(Independent_Train_dict)
Independent_t=vec.transform(Independent_Train_dict)
Independent_v=vec.transform(Independent_Val_dict)

#Standardization
scaler=StandardScaler()
scaler.fit(Independent_t)
Independent_train_sc=scaler.transform(Independent_t)
Independent_val_sc=scaler.transform(Independent_v)


#Model
model=KNeighborsRegressor(n_neighbors=15)
model.fit(Independent_train_sc, Dependent_Train)

#Predictions
VarZ_Predict=model.predict(Independent_val_sc)

#Error of model using RMSE
rmse=np.sqrt(((Dependent_Val - VarZ_Predict) ** 2).mean())

print("Error of model predictions with respect to real measurements", rmse)

Error of model predictions with respect to real measurements 0.26431968389448673

Model 2 This model only has variables that are explicitly characteristics of the wells

#Training and validation data
Training2=Wells_located_df.sample(frac=0.5)
Training2.dropna(inplace=True)
Validation2=Wells_located_df.drop(Training2.index)
Validation2.dropna(inplace=True)

variables2=["Well_Depth_Meters", "WellUse", "ActiveDays"]
Independent_Train_dict2=Training2[variables2].to_dict(orient="records")
Independent_Val_dict2=Validation2[variables2].to_dict(orient="records")

Dependent_Train2=Training2["VarZ"]
Dependent_Val2=Validation2["VarZ"]

#Dummy encoding
vec = DictVectorizer(sparse=False)
vec.fit(Independent_Train_dict2)
Independent_t2=vec.transform(Independent_Train_dict2)
Independent_v2=vec.transform(Independent_Val_dict2)

#Standardization
scaler=StandardScaler()
scaler.fit(Independent_t2)
Independent_train_sc2=scaler.transform(Independent_t2)
Independent_val_sc2=scaler.transform(Independent_v2)

#Model
model=KNeighborsRegressor(n_neighbors=15)
model.fit(Independent_train_sc2, Dependent_Train2)

#Predictions
VarZ_Predict2=model.predict(Independent_val_sc2)

#Error of model using RMSE
rmse2=np.sqrt(((Dependent_Val2 - VarZ_Predict2) ** 2).mean())

print("Error of model predictions with respect to real measurements", rmse2)

Error of model predictions with respect to real measurements 0.2323157076126313

Model 3 This model only has variables that are explicitly characteristics of the study area

#Training and validation data
Training3=Wells_located_df.sample(frac=0.5)
Training3.dropna(inplace=True)
Validation3=Wells_located_df.drop(Training3.index)
Validation3.dropna(inplace=True)

variables3=["Water_depth", "AQUIFER_NA"]
Independent_Train_dict3=Training3[variables3].to_dict(orient="records")
Independent_Val_dict3=Validation3[variables3].to_dict(orient="records")

Dependent_Train3=Training3["VarZ"]
Dependent_Val3=Validation3["VarZ"]

#Dummy encoding
vec = DictVectorizer(sparse=False)
vec.fit(Independent_Train_dict3)
Independent_t3=vec.transform(Independent_Train_dict3)
Independent_v3=vec.transform(Independent_Val_dict3)

#Standardization
scaler=StandardScaler()
scaler.fit(Independent_t3)
Independent_train_sc3=scaler.transform(Independent_t3)
Independent_val_sc3=scaler.transform(Independent_v3)

#Model
model=KNeighborsRegressor(n_neighbors=15)
model.fit(Independent_train_sc3, Dependent_Train3)

#Predictions
VarZ_Predict3=model.predict(Independent_val_sc3)

#Error of model using RMSE
rmse3=np.sqrt(((Dependent_Val3 - VarZ_Predict3) ** 2).mean())

print("Error of model predictions with respect to real measurements", rmse3)

Error of model predictions with respect to real measurements 0.23149257882329943

Model 4 This model only uses well use.

#Training and validation data
Training4=Wells_located_df.sample(frac=0.5)
Training4.dropna(inplace=True)
Validation4=Wells_located_df.drop(Training4.index)
Validation4.dropna(inplace=True)

variables4=["WellUse"]
Independent_Train_dict4=Training4[variables4].to_dict(orient="records")
Independent_Val_dict4=Validation4[variables4].to_dict(orient="records")

Dependent_Train4=Training4["VarZ"]
Dependent_Val4=Validation4["VarZ"]

#Dummy encoding
vec = DictVectorizer(sparse=False)
vec.fit(Independent_Train_dict4)
Independent_t4=vec.transform(Independent_Train_dict4)
Independent_v4=vec.transform(Independent_Val_dict4)

#Standardization
scaler=StandardScaler()
scaler.fit(Independent_t4)
Independent_train_sc4=scaler.transform(Independent_t4)
Independent_val_sc4=scaler.transform(Independent_v4)

#Model
model=KNeighborsRegressor(n_neighbors=15)
model.fit(Independent_train_sc4, Dependent_Train4)

#Predictions
VarZ_Predict4=model.predict(Independent_val_sc4)

#Error of model using RMSE
rmse4=np.sqrt(((Dependent_Val4 - VarZ_Predict4) ** 2).mean())

print("Error of model predictions with respect to real measurements", rmse4)

Error of model predictions with respect to real measurements 0.2438688822867662

Conclusions

1. The location of clusters of wells effectively occurs in areas subsiding or uplifting locally?

There are more wells in the northern area and there is more subsidence in the northern area, then considering the complete study region, YES. Considering local scales, sadly no possible to say from this analysis.

1. Are deeper wells causing more or subsidence than wells at shallow depths?

There is a really weak relation between elevation changes and well depths, yet it is positive. One well may affect areas nearby and not just under it, then I will need to see if elevation change close to each well varies more than areas without wells. Definetily there is a better relation between subsidence and well features such as well-depth, well-use and active days.

1. Is there any particular well use that may be causing more vertical changes?

Seems that there is not a strong relation between wells use and subsidence.

Challenges for future

• Is K-neighbors and appropiate model for this problem?
• I used different number of neighbors and seems that k=15 was the best, but it would be good to do a full analisys to find the optimal k (in case k-neighbors is a good model for this)